/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
	This file is part of foam-extend.

	foam-extend is free software: you can redistribute it and/or modify it
	under the terms of the GNU General Public License as published by the
	Free Software Foundation, either version 3 of the License, or (at your
	option) any later version.

	foam-extend is distributed in the hope that it will be useful, but
	WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
	General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "triSurface.H"
#include "demandDrivenData.H"
#include "IFstream.H"
#include "OFstream.H"
#include "objectRegistry.H"
#include "foamTime.H"
#include "boundBox.H"
#include "SortableList.H"
#include "PackedBoolList.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

defineTypeNameAndDebug(Foam::triSurface, 0);


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

Foam::fileName Foam::triSurface::triSurfInstance(const Time& d)
{
	fileName foamName(d.caseName() + ".ftr");

	// Search back through the time directories list to find the time
	// closest to and lower than current time

	instantList ts = d.times();
	label i;

	for (i = ts.size() - 1; i >= 0; i--)
	{
		if (ts[i].value() <= d.timeOutputValue())
		{
			break;
		}
	}

	// Noting that the current directory has already been searched
	// for mesh data, start searching from the previously stored time directory

	if (i >= 0)
	{
		for (label j = i; j >= 0; j--)
		{
			if (isFile(d.path()/ts[j].name()/typeName/foamName))
			{
				if (debug)
				{
					Pout<< " triSurface::triSurfInstance(const Time& d)"
					    << "reading " << foamName
					    << " from " << ts[j].name()/typeName
					    << endl;
				}

				return ts[j].name();
			}
		}
	}

	if (debug)
	{
		Pout<< " triSurface::triSurfInstance(const Time& d)"
			<< "reading " << foamName
			<< " from constant/" << endl;
	}

	return d.constant();
}


Foam::List<Foam::labelledTri> Foam::triSurface::convertToTri
(
	const faceList& faces,
	const label defaultRegion
)
{
	List<labelledTri> triFaces(faces.size());

	forAll (triFaces, faceI)
	{
		const face& f = faces[faceI];

		if (f.size() != 3)
		{
			FatalErrorIn
			(
				"triSurface::convertToTri"
				"(const faceList&, const label)"
			)   << "Face at position " << faceI
				<< " does not have three vertices:" << f
				<< abort(FatalError);
		}

		labelledTri& tri = triFaces[faceI];

		tri[0] = f[0];
		tri[1] = f[1];
		tri[2] = f[2];
		tri.region() = defaultRegion;
	}

	return triFaces;
}


Foam::List<Foam::labelledTri> Foam::triSurface::convertToTri
(
	const triFaceList& faces,
	const label defaultRegion
)
{
	List<labelledTri> triFaces(faces.size());

	forAll (triFaces, faceI)
	{
		const triFace& f = faces[faceI];

		labelledTri& tri = triFaces[faceI];

		tri[0] = f[0];
		tri[1] = f[1];
		tri[2] = f[2];
		tri.region() = defaultRegion;
	}

	return triFaces;
}


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::triSurface::printTriangle
(
	Ostream& os,
	const string& pre,
	const labelledTri& f,
	const pointField& points
)
{
	os  << pre.c_str() << "vertex numbers:"
		<< f[0] << ' ' << f[1] << ' ' << f[2] << endl
		<< pre.c_str() << "vertex coords :"
		<< points[f[0]] << ' ' << points[f[1]] << ' ' << points[f[2]]
		<< pre.c_str() << "region        :" << f.region() << endl
		<< endl;
}


Foam::string Foam::triSurface::getLineNoComment(IFstream& is)
{
	string line;
	do
	{
		is.getLine(line);
	}
	while ((line.empty() || line[0] == '#') && is.good());

	return line;
}


// Remove non-triangles, double triangles.
void Foam::triSurface::checkTriangles(const bool verbose)
{
	// Simple check on indices ok.
	const label maxPointI = points().size() - 1;

	forAll (*this, faceI)
	{
		const labelledTri& f = (*this)[faceI];

		if
		(
			(f[0] < 0) || (f[0] > maxPointI)
		 || (f[1] < 0) || (f[1] > maxPointI)
		 || (f[2] < 0) || (f[2] > maxPointI)
		)
		{
			FatalErrorIn("triSurface::checkTriangles(bool)")
				<< "triangle " << f
				<< " uses point indices outside point range 0.."
				<< maxPointI
				<< exit(FatalError);
		}
	}

	// Two phase process
	//   1. mark invalid faces
	//   2. pack
	// Done to keep numbering constant in phase 1

	// List of valid triangles
	boolList valid(size(), true);
	bool hasInvalid = false;

	const labelListList& fFaces = faceFaces();

	forAll (*this, faceI)
	{
		const labelledTri& f = (*this)[faceI];

		if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
		{
			// 'degenerate' triangle check
			valid[faceI] = false;
			hasInvalid = true;

			if (verbose)
			{
				WarningIn
				(
					"triSurface::checkTriangles(bool verbose)"
				)   << "triangle " << faceI
					<< " does not have three unique vertices:\n";
				printTriangle(Warning, "    ", f, points());
			}
		}
		else
		{
			// Duplicate triangle check
			const labelList& neighbours = fFaces[faceI];

			// Check if faceNeighbours use same points as this face.
			// Note: discards normal information - sides of baffle are merged.
			forAll (neighbours, neighbourI)
			{
				if (neighbours[neighbourI] <= faceI)
				{
					// lower numbered faces already checked
					continue;
				}

				const labelledTri& n = (*this)[neighbours[neighbourI]];

				if
				(
					((f[0] == n[0]) || (f[0] == n[1]) || (f[0] == n[2]))
				 && ((f[1] == n[0]) || (f[1] == n[1]) || (f[1] == n[2]))
				 && ((f[2] == n[0]) || (f[2] == n[1]) || (f[2] == n[2]))
				)
				{
					valid[faceI] = false;
					hasInvalid = true;

					if (verbose)
					{
					    WarningIn
					    (
					        "triSurface::checkTriangles(bool verbose)"
					    )   << "triangles share the same vertices:\n"
					        << "    face 1 :" << faceI << endl;
					    printTriangle(Warning, "    ", f, points());

					    Warning
					        << endl
					        << "    face 2 :"
					        << neighbours[neighbourI] << endl;
					    printTriangle(Warning, "    ", n, points());
					}

					break;
				}
			}
		}
	}

	if (hasInvalid)
	{
		// Pack
		label newFaceI = 0;
		forAll (*this, faceI)
		{
			if (valid[faceI])
			{
				const labelledTri& f = (*this)[faceI];
				(*this)[newFaceI++] = f;
			}
		}

		if (verbose)
		{
			WarningIn
			(
				"triSurface::checkTriangles(bool verbose)"
			)   << "Removing " << size() - newFaceI
				<< " illegal faces." << endl;
		}
		(*this).setSize(newFaceI);

		// Topology can change because of renumbering
		clearOut();
	}
}


// Check/fix edges with more than two triangles
void Foam::triSurface::checkEdges(const bool verbose)
{
	const labelListList& eFaces = edgeFaces();

	forAll (eFaces, edgeI)
	{
		const labelList& myFaces = eFaces[edgeI];

		if (myFaces.empty())
		{
			FatalErrorIn("triSurface::checkEdges(bool verbose)")
				<< "Edge " << edgeI << " with vertices " << edges()[edgeI]
				<< " has no edgeFaces"
				<< exit(FatalError);
		}
		else if (myFaces.size() > 2)
		{
			WarningIn
			(
				"triSurface::checkEdges(bool verbose)"
			)   << "Edge " << edgeI << " with vertices " << edges()[edgeI]
				<< " has more than 2 faces connected to it : " << myFaces
				<< endl;
		}
	}
}


// Read triangles, points from Istream
bool Foam::triSurface::read(Istream& is)
{
	is  >> patches_ >> storedPoints() >> storedFaces();

	return true;
}


// Read from file in given format
bool Foam::triSurface::read
(
	const fileName& name,
	const word& ext,
	const bool check
)
{
	if (check && !exists(name))
	{
		FatalErrorIn
		(
			"triSurface::read(const fileName&, const word&, const bool)"
		)   << "Cannnot read " << name << exit(FatalError);
	}

	if (ext == "gz")
	{
		fileName unzipName = name.lessExt();

		// Do not check for existence. Let IFstream do the unzipping.
		return read(unzipName, unzipName.ext(), false);
	}
	else if (ext == "ftr")
	{
		return read(IFstream(name)());
	}
	else if (ext == "stl")
	{
		return readSTL(name);
	}
	else if (ext == "stlb")
	{
		return readSTL(name);
	}
	else if (ext == "gts")
	{
		return readGTS(name);
	}
	else if (ext == "obj")
	{
		return readOBJ(name);
	}
	else if (ext == "off")
	{
		return readOFF(name);
	}
	else if (ext == "tri")
	{
		return readTRI(name);
	}
	else if (ext == "ac")
	{
		return readAC(name);
	}
	else if (ext == "nas")
	{
		return readNAS(name);
	}
	else
	{
		FatalErrorIn
		(
			"triSurface::read(const fileName&, const word&)"
		)   << "unknown file extension " << ext
			<< ". Supported extensions are '.ftr', '.stl', '.stlb', '.gts'"
			<< ", '.obj', '.ac', '.off', '.nas' and '.tri'"
			<< exit(FatalError);

		return false;
	}
}


// Write to file in given format
void Foam::triSurface::write
(
	const fileName& name,
	const word& ext,
	const bool sort
) const
{
	if (ext == "ftr")
	{
		return write(OFstream(name)());
	}
	else if (ext == "stl")
	{
		return writeSTLASCII(OFstream(name)());
	}
	else if (ext == "stlb")
	{
		ofstream outFile(name.c_str(), std::ios::binary);

		writeSTLBINARY(outFile);
	}
	else if (ext == "gts")
	{
		return writeGTS(sort, OFstream(name)());
	}
	else if (ext == "obj")
	{
		writeOBJ(sort, OFstream(name)());
	}
	else if (ext == "off")
	{
		writeOFF(sort, OFstream(name)());
	}
	else if (ext == "vtk")
	{
		writeVTK(sort, OFstream(name)());
	}
	else if (ext == "tri")
	{
		writeTRI(sort, OFstream(name)());
	}
	else if (ext == "dx")
	{
		writeDX(sort, OFstream(name)());
	}
	else if (ext == "ac")
	{
		writeAC(OFstream(name)());
	}
	else if (ext == "smesh")
	{
		writeSMESH(sort, OFstream(name)());
	}
	else
	{
		FatalErrorIn
		(
			"triSurface::write(const fileName&, const word&, const bool)"
		)   << "unknown file extension " << ext
			<< ". Supported extensions are '.ftr', '.stl', '.stlb', "
			<< "'.gts', '.obj', '.vtk'"
			<< ", '.off', '.dx', '.smesh', '.ac' and '.tri'"
			<< exit(FatalError);
	}
}


// Returns patch info. Sets faceMap to the indexing according to patch
// numbers. Patch numbers start at 0.
Foam::surfacePatchList Foam::triSurface::calcPatches(labelList& faceMap) const
{
	// Sort according to region numbers of labelledTri
	SortableList<label> sortedRegion(size());

	forAll (sortedRegion, faceI)
	{
		sortedRegion[faceI] = operator[](faceI).region();
	}
	sortedRegion.sort();

	faceMap = sortedRegion.indices();

	// Extend regions
	label maxRegion = patches_.size()-1;    // for non-compacted regions

	if (faceMap.size())
	{
		maxRegion = max
		(
			maxRegion,
			operator[](faceMap[faceMap.size() - 1]).region()
		);
	}

	// Get new region list
	surfacePatchList newPatches(maxRegion + 1);

	// Fill patch sizes
	forAll (*this, faceI)
	{
		label region = operator[](faceI).region();

		newPatches[region].size()++;
	}


	// Fill rest of patch info

	label startFaceI = 0;
	forAll (newPatches, newPatchI)
	{
		surfacePatch& newPatch = newPatches[newPatchI];

		newPatch.index() = newPatchI;

		label oldPatchI = newPatchI;

		// start of patch
		newPatch.start() = startFaceI;


		// Take over any information from existing patches
		if
		(
			(oldPatchI < patches_.size())
		 && (patches_[oldPatchI].name() != "")
		)
		{
			newPatch.name() = patches_[oldPatchI].name();
		}
		else
		{
			newPatch.name() = word("patch") + name(newPatchI);
		}

		if
		(
			(oldPatchI < patches_.size())
		 && (patches_[oldPatchI].geometricType() != "")
		)
		{
			newPatch.geometricType() = patches_[oldPatchI].geometricType();
		}
		else
		{
			newPatch.geometricType() = "empty";
		}

		startFaceI += newPatch.size();
	}

	return newPatches;
}


void Foam::triSurface::setDefaultPatches()
{
	labelList faceMap;

	// Get names, types and sizes
	surfacePatchList newPatches(calcPatches(faceMap));

	// Take over names and types (but not size)
	patches_.setSize(newPatches.size());

	forAll (newPatches, patchI)
	{
		patches_[patchI].index() = patchI;
		patches_[patchI].name() = newPatches[patchI].name();
		patches_[patchI].geometricType() = newPatches[patchI].geometricType();
	}
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::triSurface::triSurface()
:
	ParentType(List<Face>(), pointField()),
	patches_(0),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{}



Foam::triSurface::triSurface
(
	const List<labelledTri>& triangles,
	const geometricSurfacePatchList& patches,
	const pointField& points
)
:
	ParentType(triangles, points),
	patches_(patches),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{}


Foam::triSurface::triSurface
(
	List<labelledTri>& triangles,
	const geometricSurfacePatchList& patches,
	pointField& points,
	const bool reUse
)
:
	ParentType(triangles, points, reUse),
	patches_(patches),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{}


Foam::triSurface::triSurface
(
	const List<labelledTri>& triangles,
	const pointField& points
)
:
	ParentType(triangles, points),
	patches_(),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{
	setDefaultPatches();
}


Foam::triSurface::triSurface
(
	const triFaceList& triangles,
	const pointField& points
)
:
	ParentType(convertToTri(triangles, 0), points),
	patches_(),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{
	setDefaultPatches();
}


Foam::triSurface::triSurface(const fileName& name)
:
	ParentType(List<Face>(), pointField()),
	patches_(),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{
	word ext = name.ext();

	read(name, ext);

	setDefaultPatches();
}


Foam::triSurface::triSurface(Istream& is)
:
	ParentType(List<Face>(), pointField()),
	patches_(),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{
	read(is);

	setDefaultPatches();
}


Foam::triSurface::triSurface(const Time& d)
:
	ParentType(List<Face>(), pointField()),
	patches_(),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{
	fileName foamFile(d.caseName() + ".ftr");

	fileName foamPath(d.path()/triSurfInstance(d)/typeName/foamFile);

	IFstream foamStream(foamPath);

	read(foamStream);

	setDefaultPatches();
}


Foam::triSurface::triSurface(const triSurface& ts)
:
	ParentType(ts, ts.points()),
	patches_(ts.patches()),
	sortedEdgeFacesPtr_(nullptr),
	edgeOwnerPtr_(nullptr)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::triSurface::~triSurface()
{
	clearOut();
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::triSurface::clearTopology()
{
	ParentType::clearTopology();
	deleteDemandDrivenData(sortedEdgeFacesPtr_);
	deleteDemandDrivenData(edgeOwnerPtr_);
}


void Foam::triSurface::clearPatchMeshAddr()
{
	ParentType::clearPatchMeshAddr();
}


void Foam::triSurface::clearOut()
{
	ParentType::clearOut();

	clearTopology();
	clearPatchMeshAddr();
}


const Foam::labelListList& Foam::triSurface::sortedEdgeFaces() const
{
	if (!sortedEdgeFacesPtr_)
	{
		calcSortedEdgeFaces();
	}

	return *sortedEdgeFacesPtr_;
}


const Foam::labelList& Foam::triSurface::edgeOwner() const
{
	if (!edgeOwnerPtr_)
	{
		calcEdgeOwner();
	}

	return *edgeOwnerPtr_;
}


//- Move points
void Foam::triSurface::movePoints(const pointField& newPoints)
{
	// Remove all geometry dependent data
	deleteDemandDrivenData(sortedEdgeFacesPtr_);

	// Adapt for new point position
	ParentType::movePoints(newPoints);

	// Copy new points
	storedPoints() = newPoints;
}


// scale points
void Foam::triSurface::scalePoints(const scalar& scaleFactor)
{
	// avoid bad scaling
	if (scaleFactor > 0 && scaleFactor != 1.0)
	{
		// Remove all geometry dependent data
		clearTopology();

		// Adapt for new point position
		ParentType::movePoints(pointField());

		storedPoints() *= scaleFactor;
	}
}


// Remove non-triangles, double triangles.
void Foam::triSurface::cleanup(const bool verbose)
{
	// Merge points (already done for STL, TRI)
	stitchTriangles(points(), SMALL, verbose);

	// Merging points might have changed geometric factors
	clearOut();

	checkTriangles(verbose);

	checkEdges(verbose);
}


// Finds area, starting at faceI, delimited by borderEdge. Marks all visited
// faces (from face-edge-face walk) with currentZone.
void Foam::triSurface::markZone
(
	const boolList& borderEdge,
	const label faceI,
	const label currentZone,
	labelList& faceZone
) const
{
	// List of faces whose faceZone has been set.
	labelList changedFaces(1, faceI);

	while(true)
	{
		// Pick up neighbours of changedFaces
		dynamicLabelList newChangedFaces(2*changedFaces.size());

		forAll (changedFaces, i)
		{
			label faceI = changedFaces[i];

			const labelList& fEdges = faceEdges()[faceI];

			forAll (fEdges, i)
			{
				label edgeI = fEdges[i];

				if (!borderEdge[edgeI])
				{
					const labelList& eFaces = edgeFaces()[edgeI];

					forAll (eFaces, j)
					{
					    label nbrFaceI = eFaces[j];

					    if (faceZone[nbrFaceI] == -1)
					    {
					        faceZone[nbrFaceI] = currentZone;
					        newChangedFaces.append(nbrFaceI);
					    }
					    else if (faceZone[nbrFaceI] != currentZone)
					    {
					        FatalErrorIn
					        (
					            "triSurface::markZone(const boolList&,"
					            "const label, const label, labelList&) const"
					        )
					            << "Zones " << faceZone[nbrFaceI]
					            << " at face " << nbrFaceI
					            << " connects to zone " << currentZone
					            << " at face " << faceI
					            << abort(FatalError);
					    }
					}
				}
			}
		}

		if (newChangedFaces.empty())
		{
			break;
		}

		changedFaces.transfer(newChangedFaces);
	}
}


// Finds areas delimited by borderEdge (or 'real' edges).
// Fills faceZone accordingly
Foam::label Foam::triSurface::markZones
(
	const boolList& borderEdge,
	labelList& faceZone
) const
{
	faceZone.setSize(size());
	faceZone = -1;

	if (borderEdge.size() != nEdges())
	{
		FatalErrorIn
		(
			"triSurface::markZones"
			"(const boolList&, labelList&)"
		)
			<< "borderEdge boolList not same size as number of edges" << endl
			<< "borderEdge:" << borderEdge.size() << endl
			<< "nEdges    :" << nEdges()
			<< exit(FatalError);
	}

	label zoneI = 0;

	for (label startFaceI = 0;; zoneI++)
	{
		// Find first non-coloured face
		for (; startFaceI < size(); startFaceI++)
		{
			if (faceZone[startFaceI] == -1)
			{
				break;
			}
		}

		if (startFaceI >= size())
		{
			break;
		}

		faceZone[startFaceI] = zoneI;

		markZone(borderEdge, startFaceI, zoneI, faceZone);
	}

	return zoneI;
}


void Foam::triSurface::subsetMeshMap
(
	const boolList& include,
	labelList& pointMap,
	labelList& faceMap
) const
{
	const List<labelledTri>& locFaces = localFaces();

	label faceI = 0;
	label pointI = 0;

	faceMap.setSize(locFaces.size());

	pointMap.setSize(nPoints());

	boolList pointHad(nPoints(), false);

	forAll (include, oldFacei)
	{
		if (include[oldFacei])
		{
			// Store new faces compact
			faceMap[faceI++] = oldFacei;

			// Renumber labels for triangle
			const labelledTri& tri = locFaces[oldFacei];

			label a = tri[0];
			if (!pointHad[a])
			{
				pointHad[a] = true;
				pointMap[pointI++] = a;
			}

			label b = tri[1];
			if (!pointHad[b])
			{
				pointHad[b] = true;
				pointMap[pointI++] = b;
			}

			label c = tri[2];
			if (!pointHad[c])
			{
				pointHad[c] = true;
				pointMap[pointI++] = c;
			}
		}
	}

	// Trim
	faceMap.setSize(faceI);

	pointMap.setSize(pointI);
}


Foam::triSurface Foam::triSurface::subsetMesh
(
	const boolList& include,
	labelList& pointMap,
	labelList& faceMap
) const
{
	const pointField& locPoints = localPoints();
	const List<labelledTri>& locFaces = localFaces();

	// Fill pointMap, faceMap
	subsetMeshMap(include, pointMap, faceMap);


	// Create compact coordinate list and forward mapping array
	pointField newPoints(pointMap.size());
	labelList oldToNew(locPoints.size());
	forAll (pointMap, pointi)
	{
		newPoints[pointi] = locPoints[pointMap[pointi]];
		oldToNew[pointMap[pointi]] = pointi;
	}

	// Renumber triangle node labels and compact
	List<labelledTri> newTriangles(faceMap.size());

	forAll (faceMap, facei)
	{
		// Get old vertex labels
		const labelledTri& tri = locFaces[faceMap[facei]];

		newTriangles[facei][0] = oldToNew[tri[0]];
		newTriangles[facei][1] = oldToNew[tri[1]];
		newTriangles[facei][2] = oldToNew[tri[2]];
		newTriangles[facei].region() = tri.region();
	}

	// Construct subsurface
	return triSurface(newTriangles, patches(), newPoints, true);
}


void Foam::triSurface::write
(
	const fileName& name,
	const bool sortByRegion
) const
{
	write(name, name.ext(), sortByRegion);
}


void Foam::triSurface::write(Ostream& os) const
{
	os  << patches() << endl;

	//Note: Write with global point numbering
	os  << points() << nl
		<< static_cast<const List<labelledTri>&>(*this) << endl;

	// Check state of Ostream
	os.check("triSurface::write(Ostream&)");
}


void Foam::triSurface::write(const Time& d) const
{
	fileName foamFile(d.caseName() + ".ftr");

	fileName foamPath(d.path()/triSurfInstance(d)/typeName/foamFile);

	OFstream foamStream(foamPath);

	write(foamStream);
}


void Foam::triSurface::writeStats(Ostream& os) const
{
	// Unfortunately nPoints constructs meshPoints() so do compact version
	// ourselves.
	PackedBoolList pointIsUsed(points().size());

	label nPoints = 0;
	boundBox bb = boundBox::invertedBox;

	forAll (*this, triI)
	{
		const labelledTri& f = operator[](triI);

		forAll (f, fp)
		{
			label pointI = f[fp];
			if (pointIsUsed.set(pointI, 1))
			{
				bb.min() = ::Foam::min(bb.min(), points()[pointI]);
				bb.max() = ::Foam::max(bb.max(), points()[pointI]);
				nPoints++;
			}
		}
	}

	os  << "Triangles    : " << size() << endl
		<< "Vertices     : " << nPoints << endl
		<< "Bounding Box : " << bb << endl;
}


// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //

void Foam::triSurface::operator=(const triSurface& ts)
{
	List<labelledTri>::operator=(ts);
	clearOut();
	storedPoints() = ts.points();
	patches_ = ts.patches();
}


// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //

Foam::Ostream& Foam::operator<<(Ostream& os, const triSurface& sm)
{
	sm.write(os);
	return os;
}


// ************************************************************************* //
